Numerical simulations of directed sandpile models 
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Directed sandpile models with different toppling rules are studied by means of numerical sim- 
ulations in two dimensions, with the purpose of determining the different universality classes. It 
is concluded that the random-threshold directed model is in the same universality of the Manna 
directed model where multiple toppling events plays a determinant role. The BTW model with 
uniform and noisy driving are found to display the same critical behavior. Moreover it is observed 
that the Zhang model does not satisfy simple finite size scaling. 
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I. INTRODUCTION 

The classification of sandpile models in their different 
universality class has been a topic of intensive research 
in the last decade pL However, most of the work has 
been devoted to models with undirected toppling while 
their corresponding directed variants has been less stud- 
ied. From the theoretical side we count with the exact so- 
lution for the directed BTW model ^ obtained by Dhar 
and Ramaswamy ||, which can be taken as a reference 
for numerical simulations. On the other hand, Pastor- 
Satorras and Vespignani has recently reported nu- 
merical simulations for the Manna model || , which puts 
clear evidence about the classification of the BTW and 
Manna directed models in different universality classes. 
Other studies include a directed model with a probabilis- 
tic toppling ||] and a recent report concerning the effect 
of local dissipation on the SOC state 0. 

The analysis of other directed models, the non Abelian 
Zhang model for instance, is of great interested and may 
give light to the study of their corresponding undirected 
variants. With this scope three different directed models 
are investigated by means of numerical simulations. This 
includes the well known Zhang model Q, the random 
threshold model (RT) Jj| , and the BTW model under an 
uniform driving. The influence of a uniform driving in the 
Zhang model has been already discussed in the literature 
although some aspects are still not clear. However, 
the same analysis has not been made for models with a 
discrete toppling rule, like the BTW model, where the 
energy transfer always takes place in discrete units. In 
this direction Narayan and Middleton Jll]] suggested that 
the BTW model under noisy and uniform driving has the 
same critical behavior. 

From the numerical data it is concluded that the 
Manna and RT directed models belong to the same uni- 
versality class, in agreement with the gen eral believe for 
the corresponding undirected variants |l^Jl^] . Moreover, 
the data obtained for the Zhang model does not satisfy 
finite-size scaling which suggest that this model is in dif- 



ferent universality class from that of the models men- 
tioned above. 

In the case of the BTW model under uniform driving it 
is shown, after some algebra with the toppling operators 
||, that its evolution is periodic in time with a period 
scaling linearly with system size. In spite of this period- 
icity, which is not present in the original model with noisy 
driving || , the statistics of the avalanches is found to be 
practically identical to its noisy driving counterpart. 



II. MODELS AND SIMULATIONS 
A. Models with noisy driving 

Consider a square lattice of L 2 sites labeled by index 
(h j) (h j = 1 , - - - , N) and assign a variable Zij to each 
of them, be continuous or discrete and may have 

different interpretations depending on the system one is 
modeling. It will be referred here as the energy storage 
at the corresponding site. The geometry used is shown 
in fig. [l], in which a site can transfer energy only to its 
three downward nearest neighbors (nn). In the horizon- 
tal direction periodic boundary conditions are considered 
while the downward boundary is taken open. This ge- 
ometry allows a natural implementation of the Manna 
toppling rule Q . 

One motivation for the use of this geometry was given 
in it is just introduced to allow the implementation 
of the Manna toppling rule. 

To completely define a model one should specify the 
initial condition and the evolution rules (addition of en- 
ergy and toppling) of the sandpile cellular automaton 
model. A threshold z c is considered in such a way that 
sites with z < z c are say to be stable and their energy 
remains constant, while those with z > z c are say to be 
active and topple transferring energy to their downward 
nn. First the usual noisy addition of energy is considered. 
In this case, if all sites are stable a unit of energy is added 
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to a site selected at random. Then the system is updated 
in parallel using the toppling rule until all sites are stable. 
The number of toppling events required to drive the sys- 
tem to an stable configuration is the size of the avalanche 
and is denoted by s. On the other hand, the number of 
steps (parallel updates) required is its duration and is 
denoted by T. Since the driving acts at random after 
some avalanches the system "forget" its initial condition 
and reaches a stationary state. In other words the initial 
condition is irrelevant. 

Models will differ from each other depending on the 
specific toppling rule one implements. Here the following 
toppling rules are considered, 

• BTW: z c = 3, Zij —> z i3 - 3 and Zkj-i -> zjy-i + 1 
(k = i — + 1); 

• Manna: z c — 2, z^ — ► Zy — 2 and Zkj-i — » 
%kj—i + <5fc (k = i — + where 5k can take the 
values 0, 1, 2 at random but with the constraint of 
conservation ^2 k 5k = 2; 

• Zhang: z c = 1, z {j — > and Zkj-i —> Zjy_i + zy/3 
(k = i — 1,m + 1); 

• RT: z c takes the values 3 and 4 at random after 
each toppling, z^ — > Zjj — 3 and Zkj~i — ► Zfcj_i + 1 
(fc = i- 1,« + 1). 

B. BTW model under uniform driving 

In the noisy driving described above the addition of 
energy takes place at one site selected at random. How- 
ever, there are many situations where a uniform driving 
in which the energy at all sites increases in the same 
amount becomes more realistic. Examples can be found 
in earthquake dynamics [fiol , interface depinning | fj"i"| and 
also in some experimental setups to for granular mate- 
rials fyj. This type of driving has been investigated in 
models with a toppling rule similar to that of the Zhang 
models but with local dissipation fljof . 

In the case of the BTW model we should be careful 
when introducing a global driving. If as usual Zij is an 
integer variable and the energy at all sites is increased at 
a constant rate c then many sizes will reach the threshold 
energy at the same time and, therefore, many avalanches 
will start at different points of the lattice leading to the 
superposition of avalanches. 

This problem can be solved considering a continuum 
energy profile. This is still not enough because if all sites 
start with a discrete energy it will remain discrete for- 
ever. We are thus forced to consider a continuum initial 
profile Zij(0). Then as it was already shown by Narayan 
and Middleton flTf the continuum addition of energy can 
be replaced by a sequential addition of energy. For sim- 
plicity consider the low disorder regime where zij (0) < 1 
for all sites. For the analysis below is irrelevant if the 
initial energy profile is displaced uniformly at all sites. 



Now, suppose that the energy increases at rate c at all 
sites. An example is shown in fig. || for a lattice made of a 
horizontal line of three sites. Notice that in this case one 
has only input of energy coming form the driving field 
and output dissipation under toppling, a simplification 
considered for illustrative purposes. 

In the continuum time scale the energy increases lin- 
early until it reaches the threshold where the site topples 
and its energy decreases by 3 (BTW toppling rule). But 
the system can be also monitored in a discrete time and 
energy scales. In these scales, at step t = all sites 
have energy 0. In step 1, 2 and 3 sites 1, 2 and 3 receives 
one unit of energy. Then in subsequent steps the same se- 
quence of addition is repeated. The order in the sequence 
of addition is clearly determined by the initial condition 
and all sites should receive a unit of energy before the 
first site of the second receives the second grain. 

Now consider an square lattice, where sites can also 
receive energy from nearest neighbors in the layer above. 
The picture will not change in relation to the addition 
of energy from the external field. In the BTW model 
the energy is transferred in discrete units and, therefore, 
the toppling only modifies the integer part of z with no 
modification of the sequence of addition. This is a fun- 
damental difference with the Zhang toppling rule which 
not only involve the integer par of the energy but on top- 
pling all the energy at the active site is transferred. The 
consequences derived from this periodic sequential driv- 
ing is investigated below, using the formalism introduced 
by Dhar et al ||. 

Let be dij the operator which add a particle at site 
and lead the system relax to an equilibrium posi- 
tion ||. After N steps all sites receives one, and only 
one, unit of energy in certain order determined at t = 0. 
Thus, if at time t we have a configuration C(t) then at 
time t + N we will obtain the configuration 

L L 

C(t + JV) = JJJJoyC(t). (1) 

i=lj=l 

The order in which the string of operators appears in this 
equation is irrelevant because the operators aij commute 
among them self. 

Applying these string three times it results that 

L L 

c(i+3A0=nn4 c w- ( 2 ) 

i=U=l 

This expression can be simplified using the following 
property of the toppling operators || 

3 _ \ a i-ij-i a ij-\Q-i+ij-i, for j < L, , , 

fl y _ \ 1, for j = L. [6) 

The first equality express the fact that the addition of 
three grains at a site with j < L, makes this site 

active transferring one grain to each of its downward nn. 
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The second one applies for the boundary sites which af- 
ter receiving three grains become active dissipating these 
three grains through the boundary and, therefore, leaving 
the energy configuration invariant. 

Starting at layer j — 1 all the operators an are elim- 
inated using eq. (|3j). This will increase the power of 
operators aa in eq. (^) by 3. The same procedure is ap- 
plied to the second, third ... L — l layer finally resulting 

L 

C(t + 3N) = l[a%C(t). (4) 

i=i 

The application of the operator an, three consecutive 
times lead its energy invariant and, therefore, eq. (Q) 
is reduced to 

C(t + 3N) =C(t). (5) 

Hence the evolution of the energy profile is periodic with 
period 3N. 

This property is not observed in the noisy driving case 
where the randomness introduced by the driving field 
makes the dynamics Markovian ||. Nevertheless, as it is 
shown in the next section, the statistical properties of the 
avalanches in the BTW directed model are independent 
of the driving mechanism. 

C. Numerical simulations and discussion 

Numerical simulations of the BTW, Manna, Zhang and 
RT models with directed toppling were performed. In all 
cases lattice sizes ranging from L = 64 to L = 2048 where 
used. The numerical results obtained for the BTW and 
Manna models are taken only as a reference because we 
count with the largest scale simulations reported in 
(up to L = 6400). 

Noisy driving: Starting from an initial flat profile all 
systems were updated until they reach the stationary 
state. After that statistics over 10 8 avalanches was taken, 
recording the avalanche sizes and durations. 

Uniform driving: the evolution in time of the energy 
profile is periodic and, therefore, average was taken over 
the period 37V. Different initial conditions where sim- 
ulated using different permutations of the sequence of 
addition of energy. 

To extract the scaling exponents we use the moment 
analysis technique ]15l| . The q moment of the probability 
density p x (x) of a magnitude x is defined by 

(x q ) = Jdx Px (x)x q . (6) 

where x — s,T. As defined above s and T are the 
avalanche size and duration, i.e. the number of top- 
pling events and parallel updates, respectively, required 
to drive the system to an stable configuration. 

If the hypothesis of finite-size scaling is satisfied, that 
is the distribution of avalanche size and duration can be 



written in the form p x (x) — x~ T " f x (x/L^ :c ), then the q 
moment scales with system size according to the power 
law 

(x q ) ~ L a *( q \ (7) 

with 

<r x (q) = p x (i - t x ) + p x q, (8) 

where f} s = D and /3t = z are effective dimensions 
which characterize how the cutoffs of the distribution of 
avalanche sizes and durations, respectively, scales with 
system size. On the other hand r x is the power low expo- 
nent which can be measured in the scaling region before 
the finite-size cutoffs. 

The plot of o~ s (q) and o~t{o) vs. q is shown in figs. 
H for different directed models, and ||, respectively. If 
two models belong to the same universality class then 
the linear part of the plot should overlap. Based on this 
argument it is then concluded that the RT model be- 
long to the same universality class of the Manna model. 
A more quantitative comparison can be seen in table || 
where the exponents computed here for the RT model 
are compared with those reported in Q for the Manna 
model. The scaling exponents are found in very good 
agreement within the numerical error. 

If the hypothesis of finite size scaling is valid then one 
can take the scaling exponents obtained from the moment 
analysis and plot the different distributions on rescaled 
variables in such a way that curve for different system 
sizes overlap. This is done in figs. || and || for the RT 
model resulting in a very good data collapse, as it has 
been also observed for the directed Manna model 

On the other hand, one cannot distinguish between the 
curve for the BTW model with noisy or uniform driv- 
ing, leading to the same scaling exponents. Thus, the 
periodicity introduced by the uniform driving carry no 
consequence for the critical behavior of the BTW model. 
Hence, the noisy driving can be substituted by a uniform 
driving together with an initial random energy profile. 
This will correspond in an interface depinning descrip- 
tion, the number of toppling events playing the role of 
the interface height, to a columnar disorder. A similar 
conclusion was obtained by Lauritsen and Alava using a 
different argument p3| . 

The things becomes less clear when analyzing the 
Zhang model. In this case from the moment technique it 
results that D w 1.55, z k, 1.03, r s rj 1.31 and r t « 1.53. 
These exponents define by them self a new universality 
class. However, the moment analysis technique is based 
on the hypothesis of finite size scaling which in this case 
is not satisfied. This fact becomes clear in figs, [z] and ||, 
where the data collapse is shown, revealing that in this 
case the finite-size scaling is not satisfied. Deviations are 
observed not only for the smallest avalanches but also 
for the largest avalanches where the finite size scaling is 
expected to be better. 
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The anomalies observed for the Zhang model are as- 
sociated with the existence of huge avalanches which 
practically empties the system. After one of these huge 
avalanches the system needs some time to reach again 
the critical state. This means that the mean energy 
of the system displays strong fluctuations and, there- 
fore, the overall avalanche statistics is given by the small 
avalanches taking place during the accumulation of new 
grains and these huge avalanches. This picture is illus- 
trated in fig. ^ where the fraction of avalanches of size s 
is plotted. It is characterized by a rounded peak at the 
largest avalanche sizes which shifts with lattice size. On 
the other hand, the other part of the distribution cannot 
be fitted by a single power law. 

The classification of the Manna and RT directed mod- 
els in the same universality class is in agreement with 
a similar report for the corresponding undirected vari- 
ants |12|. Thus, there should be some common element 
in these models, which is off course not present in the 
BTW model. A clue was given in B related with the 
possibility of multiple toppling events. In this final part 
of this section we discuss this statement in more detail. 

In the directed BTW model the cluster of sites which 
topples within an avalanche is compact and these sites 
topple only once. On the contrary, Pastor-Satorras and 
Vespignani Q observed that in the directed Manna model 
the cluster of sites touched by the avalanche is still com- 
pact but each site participating in the avalanche can top- 
ple more than once. If the existence of multiple toppling 
events is the property that puts the Manna model in a 
different universality class then a similar behavior should 
be observed in the present simulations of the directed RT 
model. 

A decomposition of the sites participating in an 
avalanche based on the number of toppling events per- 
formed at these sites is shown in fig. 9, for the case of 
the directed random-threshold model. In this particular 
realization the cluster of sites touched by the avalanche 
is decomposed in three sub clusters where sites have top- 
pled one, two and three times. The fraction of sites top- 
pling three times is small but the one with two toppling is 
comparable with that of one. In general it was observed 
that in large avalanches the fraction of sites which topple 
one and two times are of the same order and, therefore, 
multiple toppling events are relevant. 

In the case of undirected models it is known that multi- 
ple toppling events are present even in the BTW model, 
which leads to the decomposition of the avalanches in 
waves [Tq| . However, their origin is different than in di- 
rected models. In the undirected BTW model a site may 
topple more than once because after a first toppling (let 
say at step t) it is possible that all its neighbors become 
active and topple (at step t + 1) and, therefore, the site 
will again be active and topple (at step t + 2). In the 
decomposition of waves one apply the toppling rule to 
all sites until they are stable before toppling the initial 
active site a second, third, ... time, generating in this 
way the first, second, ... wave. 



One may think in applying a similar approach for the 
avalanches in the Manna and RT directed models, de- 
composing the avalanche as a superposition of waves. A 
fundamental property of the waves is that within it sites 
can toppling only once, otherwise the concept is useless. 
Below its is shown that such a decomposition is not pos- 
sible in the Manna and RT directed models, a least not 
in such a simple way. 

Let us analyze in detail how a multiple toppling event 
can be generated in the RT directed model. Suppose 
the lattice has a configuration where a site has height 3 
and threshold 4 and its three upward nearest neighbors 
are active. Then in the next step the site will receive 
three grains, one per active neighbor, taking an energy 
3 + 3 = 6 > 4, becoming active. After toppling the en- 
ergy will decrease to 6 — 3 = 3 and a new threshold is 
assigned. But the new threshold can be either 3 or 4. 
If it is 4 the site will be stable but if it is 3 it will be 
still active and topple in the next step. Since in the par- 
ticular model considered here the two threshold are se- 
lected with equal probability then the multiple toppling 
can take place with the same probability than the sin- 
gle one, which explains the previous observation that in 
large avalanche the fraction of two-toppling events at the 
same site are of the order of the one-toppling one. 

During the evolution of an avalanche which started at 
layer jo is possible that a site at a layer below jj > jo 
needs two consecutive toppling to be stable. Thinking 
in a decomposition in waves one can delay the second 
toppling until all the sites below are stable (first wave) 
and then topple the site the second time generating the 
second wave. However, during the first wave is possible 
that a site at a deeper layer ji > ji > jo also needs 
two toppling to become stable and, therefore, the first 
wave has to be decomposed in sub-waves where sites 
topple only once. The same process may occur even 
at deeper layer thus generating a hierarchical structure 
of sub-waves. Hence, the decomposition of avalanche in 
waves in these models lead to a more complex structure 
which nevertheless may be exploited to obtain some es- 
timate of the scaling exponents. This is nevertheless out 
of the scope of this work. 



III. SUMMARY AND CONCLUSIONS 

Directed sandpile models with different toppling rule 
has been studied by means of numerical simulations, 
with the purpose of determining the different universal- 
ity classes. To extract the scaling exponents the moment 
analysis technique was used and the resulted exponents 
were latter corroborated by finite size scaling of the dis- 
tribution of avalanche size and duration. 

The numerical analysis reveals that the introduction 
of a uniform driving in the BTW directed model does 
not change the critical properties. The evolution in time 
of the energy profile is in this case periodic with a pe- 
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riod which scales linearly with the system size. In spite 
of this periodicity the avalanche distributions are practi- 
cally identical to that obtained for the same model but 
with the usual noisy driving. 

It is concluded that the Manna and RT models are 
in the same universality, where multiple toppling events 
appear to be a fundamental property. The existence of 
multiple toppling events leads to a decomposition of the 
avalanche in a hierarchical structure of waves which my 
be a starting point for future research. 

Finally, it is observed that the avalanches in the di- 
rected Zhang model displays a complex structure which 
does not satisfy the finite-size scaling hypothesis. It is 
given by the superposition of huge avalanches involving a 
large dissipation of energy through the boundary a small 
avalanches taking place during the accumulation of en- 
ergy. 
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FIG. 1. Geometry used in the simulations. A square lat- 
tice oriented from top to bottom with each site having three 
upward and three downward nearest neighbors. 



Continuous time 
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Discrete time 

FIG. 2. Schematic mapping of the BTW model under uni- 
form driving with a continuum energy profile to the same 
model with a sequential driving and a discrete energy profile. 
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FIG. 3. The exponents a s (q) as a function of g. From top 
to bottom appear the curve for the Manna, RT, Zhang and 
BTW directed models. In all cases except the one labelled 
by CD (continuous driving) a noisy driving was used. In the 
case of the BTW directed model one can not distinguish the 
curves for noisy and uniform or continuous driving. 
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FIG. 4. The exponents (JT{q) as a function of q. From top 
to bottom appear the curve for the Zhang, BTW, RT and 
Manna directed models. In the case of the BTW directed 
model one can not distinguish the curves for noisy and uni- 
form driving (CD), models. 
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FIG. 5. Finite size scaling of the integrated distribution of 
avalanche sizes (avalanches of size greater that s) of the RT 
directed model, using the scaling exponents reported in tab. 
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FIG. 6. Finite size scaling of the integrated distribution of 
avalanche durations (avalanches of duration larger than T) of 
the RT directed model, using the scaling exponents reported 
in tab. M. 
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FIG. 7. Finite size scaling of the distribution of avalanche 



sizes of the Zhang directed model, using D 

T 3 = 1.31. 
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FIG. 8. Finite size scaling of the distribution of avalanche 
durations of the Zhang directed model, using z — 1.03 and 
Tt = 1.53. 
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FIG. 9. Fraction of avalanches of size s in the Zhang model 
for three different lattice sizes. Notice the existence of a 
rounded peak for the largest avalanche sizes with shift with 
system size. 
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FIG. 10. Number of toppling events at each lattice size in 
an avalanche taking place in a system of size L — 64, using the 
RT toppling rule. Toppling takes place from small to large j. 
Three layers are clearly observed, corresponding to one, two 
and three toppling. 



Model 


D 


z 




Tt 


BTW 1 


3/2 


1 


4/3 


3/2 


MannaT 


| 1.75(1) 


0.99(1) 


1.43(1) 


1.74(4) 


RT 


1.73(2) 


0.99(3) 


1.44(2) 


1.70(4) 



TABLE I. The scaling exponents for the BTW, Manna and 
RT directed models. Those for the BTW model are exact 
while the others are numerical estimates. 
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